#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <iostream>
#include <cmath>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/io.hpp>

const double PI = (std::atan(1.0) * 4.0);

double source(double x, double y);

double u_b(double x, double y);

double source(double x, double y)
{
    return 2.0 * PI * PI * std::sin(PI * x) * std::sin(PI * y);
};

double u_b(double x, double y)
{
    return std::sin(PI * x) * std::sin(PI * y);
};


int main(int argc, char *argv[])
{

    using namespace boost::numeric::ublas;

    double x0 = 0.0;
    double y0 = 0.0;
    double x1 = 1.0;
    double y1 = 1.0;

    long N = 100;
    double h = (x1 - x0) / (N - 1);
    long m = N - 2;
    long n = m * m;

    compressed_matrix<double> A (n, n);

    typedef boost::numeric::ublas::compressed_matrix<double>::iterator1 it1_t;
    typedef boost::numeric::ublas::compressed_matrix<double>::const_iterator1 it1_t2;
    typedef boost::numeric::ublas::compressed_matrix<double>::iterator2 it2_t;
    typedef boost::numeric::ublas::compressed_matrix<double>::const_iterator2 it2_t2;

    vector<double> f (n);

    for (long i = 0; i < n; i++)
	f (i) = 0;

    vector<double> u (n);

    for (long i = 0; i < m; i++)
	for (long j = 0; j < m; j++)
	{
	    long idx = i * m + j;
	    double lx = (j + 1) * h;
	    double ly = (i + 1) * h;
	    A (idx, idx) = 4.0;
	    f (idx) = source(lx, ly) * h * h;
	    if (j != 0)
		A (idx, idx - 1) = -1.0;
	    else
		f (idx) += u_b(lx - h, ly);
	    if (j != m - 1)
		A (idx, idx + 1) = -1.0;
	    else
		f (idx) += u_b(lx + h, ly);
	    if (i != 0)
		A (idx, idx - m) = -1.0;
	    else
		f[idx] += u_b(lx, ly - h);
	    if (i != m - 1)
		A (idx, idx + m) =  -1.0;
	    else
		f (idx) += u_b(lx, ly + h);
	}
    
    for (it1_t t1 = A.begin1(); t1 != A.end1(); t1++)
	for (it2_t t2 = t1.begin(); t2 != t1.end(); t2++)
	    std::cout << *t2 << std::endl;
    return 0;
};
